Models with covarying parameters

Jesse Brunner

Example: Growing virus isolates

Recall from lab that exponentially growing virus, \(V\) should look like this: \[ V(t) = V_0 e^{rt} \] which, on the log scale, is a nice linear model: \[ \log[V(t)] = \log(V_0) + rt \] Our study:

  • ten strains of virus grown in cell culture
  • take samples on days 1, 3, and 5 and titer the virus
  • expect strains to vary in their growth rates
  • Moreover, we expect that some viruses might be very infectious, and so have a high \(V_0\), but this might come at the expense of a low growth rate, \(r\).

Simulating our example: how to deal with covariances

Need to simulate data where \(r\) and \(V_0\) are correlated for each strain.

r_mu <- 1.5 # mean growth rate
r_sig <- 0.5 # sd of growth rate
V0_mu <- 4.6 # mean initial virus population = log(100) 
V0_sig <- 1.5 # sd of initial virus population
rho <- -0.65 # correlation between V0 and r

# vector of means
Mu <- c(V0_mu, r_mu)
# vector of standard deviations
sigmas <- c(V0_sig, r_sig)

Option 1 to build correlation matrix:

# covariance of the two parameters
(cov_V0r <- r_sig * V0_sig * rho)
[1] -0.4875
# covariance matrix
(Sigma <- matrix( c(V0_sig^2, cov_V0r, 
                    cov_V0r, r_sig^2), 
                  ncol=2, byrow = TRUE))
        [,1]    [,2]
[1,]  2.2500 -0.4875
[2,] -0.4875  0.2500

Option 2 to build correlation matrix:

# first make correlation matrix
(Rho <- matrix( c(1, rho, 
                  rho, 1), 
                ncol=2, byrow=TRUE))
      [,1]  [,2]
[1,]  1.00 -0.65
[2,] -0.65  1.00
# Then matrix multiply to get covariance matrix
(Sigma <- diag(sigmas) %*% Rho %*% diag(sigmas))
        [,1]    [,2]
[1,]  2.2500 -0.4875
[2,] -0.4875  0.2500

\[ \begin{bmatrix} a & b \\ c & d \end{bmatrix} \times \begin{bmatrix} e & f \\ g & h \end{bmatrix} = \begin{bmatrix} ae +bg & af + bh \\ ce + dg & cf + dh \end{bmatrix} \]

Simulating our example: how to deal with covariances

Simulate parameters by group

nS <- 20 # number of strains

params <- MASS::mvrnorm(n=nS, mu=Mu, Sigma=Sigma)
colnames(params) <- c("V0", "r")
(params <- as.data.frame(params))
          V0         r
1  5.6997115 0.8928427
2  4.4961836 1.4243006
3  6.3765896 1.1834381
4  4.9737238 1.9757917
5  7.4924464 1.2853827
6  4.7185409 1.1903997
7  3.4052933 1.9689721
8  2.7186144 1.4827463
9  6.1552155 0.7785731
10 5.1505684 1.3216898
11 6.0818282 1.5490029
12 4.4921862 1.2598660
13 4.3676307 1.3278242
14 0.7001111 2.2344075
15 5.0045218 1.4107406
16 5.9830801 1.4213032
17 3.8837860 1.5354338
18 4.3212154 1.0806057
19 5.1253671 2.1950777
20 3.3916591 1.6799052
plot(r ~ V0, data = params)

Simulating our example: how to deal with covariances

Finally time to simulate observations!

time <- c(1,3,5) # time points
sigma <- 0.25 # observation error

df <- expand.grid(ID = 1:nS,
                  time = time)
df$V <- rnorm(n=nrow(df), 
                 mean=params$V0[df$ID] + params$r[df$ID]*df$time, 
                 sd=sigma)

Simulating our example: how to deal with covariances

Analysis goal

Our research questions are:

  1. Do virus strains have substantially different \(r\)’s and \(V_0\)’s
  2. Are the \(r\)’s and \(V_0\)’s negatively correlated?

How would you analyze these data?

Varying slopes and intercepts model

a & b vary by clusters, but are independent \[ \begin{align} \log(V_i) &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= a_{\text{Strain}[i]} + b_{\text{Strain}[i]]}\times \text{time} \\ a_{\text{Strain}} &\sim \text{Normal}(\mu_{a},\sigma_a) \\ b_{\text{Strain}} &\sim \text{Normal}(\mu_{b}, \sigma_b) \\ \mu_{a} &\sim \text{Normal}(4,1.5) \\ \sigma_{a} &\sim \text{Exponential}(2) \\ \mu_{b} &\sim \text{Normal}(1,1) \\ \sigma_{b} &\sim \text{Exponential}(3) \\ \sigma &\sim \text{Exponential}(2) \end{align} \]

a & b vary by clusters, but are correlated $$ \[\begin{align} \log(V_i) &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= a_{\text{Strain}[i]} + b_{\text{Strain}[i]]}\times \text{time} \\ \left[\begin{matrix} a_{\text{Strain}} \\ b_{\text{Strain}} \end{matrix} \right] &\sim \text{MVNormal}\left(\left[\begin{matrix} \mu_{a} \\ \mu_{b} \end{matrix} \right],\Sigma \right) \\ \Sigma &= \left( \begin{matrix} \sigma_a & 0 \\ 0 & \sigma_b \end{matrix} \right) \text{Rho} \left( \begin{matrix} \sigma_a & 0 \\ 0 & \sigma_b \end{matrix} \right) \\ \mu_{a} &\sim \text{Normal}(4,1.5) \\ \sigma_{a} &\sim \text{Exponential}(2) \\ \mu_{b} &\sim \text{Normal}(1,1) \\ \sigma_{b} &\sim \text{Exponential}(3) \\ \text{Rho} &\sim \text{LKJcorr}(2) \\ \sigma &\sim \text{Exponential}(2) \end{align}\] $$

Varying slopes and intercepts model, in ulam()

a & b vary by clusters, but are independent

m1 <- ulam(
  alist(
    V ~ dnorm(mu, sigma), 
    mu <- a[ID] + b[ID]*time,
    
    # priors
    a[ID] ~ dnorm(a_mu, a_sd),
    a_mu ~ dnorm(4,1.5),
    a_sd ~ dexp(2),
    
    b[ID] ~ dnorm(b_mu, b_sd),
    b_mu ~ dnorm(1,1),
    b_sd ~ dexp(3),
    
    sigma ~ dexp(2)
  ), data = df
)

a & b vary by clusters, but are correlated

m2 <- ulam(
  alist(
    V ~ dnorm(mu, sigma), 
    mu <- a[ID] + b[ID]*time,
    
    # priors
    c(a, b)[ID] ~ multi_normal( c(a_mu, b_mu), Rho, sigma_ID),
    a_mu ~ dnorm(4,1.5),
    b_mu ~ dnorm(1,1),
    sigma_ID ~ dexp(2),
    
    Rho ~ lkj_corr(2),
    sigma ~ dexp(2)
  ), data = df
)

Varying slopes and intercepts model, in ulam()

a & b vary by clusters, but are independent

precis(m1, depth=2)
           mean         sd      5.5%     94.5%    n_eff     Rhat4
a[1]  5.6197359 0.28164215 5.1356684 6.0741868 449.9927 0.9981760
a[2]  4.7118293 0.28615075 4.2662249 5.1606819 487.9817 0.9980835
a[3]  6.2525709 0.26581581 5.8259082 6.6684760 375.2530 0.9980828
a[4]  4.6659927 0.27625705 4.2391181 5.1218754 445.6447 0.9982450
a[5]  7.4975931 0.28200501 7.0029107 7.9305953 446.4993 1.0010159
a[6]  4.6477172 0.28387578 4.2254237 5.1334320 515.7942 0.9989191
a[7]  3.8892686 0.28299326 3.4693857 4.3807230 557.8814 1.0045932
a[8]  2.9734650 0.27109802 2.5508475 3.3999261 376.7245 1.0004947
a[9]  5.6376343 0.29488779 5.1884568 6.0771947 446.8889 1.0029273
a[10] 5.1894291 0.27536434 4.7141119 5.6159996 438.3649 1.0018388
a[11] 5.8634807 0.27772760 5.3923784 6.3354943 451.7656 0.9982016
a[12] 4.1641750 0.29046496 3.7086837 4.6194622 470.7299 0.9980760
a[13] 4.3872350 0.27518780 3.9770622 4.8469286 505.0801 0.9989898
a[14] 1.2771100 0.30963280 0.8441708 1.7853930 307.5220 0.9998305
a[15] 5.2451792 0.25533451 4.8367318 5.6186080 401.0076 0.9983295
a[16] 5.8120824 0.26097571 5.3944576 6.2242712 642.1812 0.9980392
a[17] 4.0266713 0.28655454 3.5669885 4.5023922 467.8292 0.9994088
a[18] 4.2720804 0.26584935 3.8388472 4.6814475 404.6871 0.9988240
a[19] 5.7211973 0.27967798 5.2838661 6.1614799 610.1397 0.9992123
a[20] 3.5095258 0.28534628 3.0684600 3.9735664 372.2920 1.0034292
a_mu  4.7252093 0.29493601 4.2705324 5.2030886 697.2526 0.9988166
a_sd  1.3396701 0.20974641 1.0563460 1.6870729 508.1255 1.0015487
b[1]  1.0294046 0.08191225 0.8973560 1.1581166 409.2617 0.9993012
b[2]  1.4568568 0.08550066 1.3160524 1.5857219 499.0913 0.9980833
b[3]  1.2596826 0.08064037 1.1295799 1.3928521 370.2664 0.9984688
b[4]  2.0746709 0.07813137 1.9473135 2.1932221 372.4381 0.9987550
b[5]  1.2281400 0.08432477 1.0976268 1.3665780 364.2279 1.0003461
b[6]  1.2463987 0.08155160 1.1141551 1.3718674 494.7723 0.9996309
b[7]  1.8302930 0.08344857 1.6939483 1.9615252 504.1925 1.0054988
b[8]  1.3686721 0.07815777 1.2420503 1.4925317 372.4234 1.0021726
b[9]  0.8990466 0.08866126 0.7627035 1.0444876 526.6238 1.0018244
b[10] 1.3045427 0.08006247 1.1831102 1.4277186 468.2089 0.9998287
b[11] 1.5981687 0.08279476 1.4689424 1.7377499 498.9307 0.9983048
b[12] 1.2996811 0.08305029 1.1676456 1.4362483 467.3709 0.9981242
b[13] 1.3551176 0.08227767 1.2229513 1.4774558 505.2874 0.9988714
b[14] 2.0515093 0.09284407 1.9039347 2.1831765 302.6048 1.0009567
b[15] 1.3380622 0.07448518 1.2197519 1.4601463 425.0079 0.9984582
b[16] 1.4086105 0.07839334 1.2891097 1.5310297 680.2872 0.9981991
b[17] 1.4665637 0.08284514 1.3389944 1.5871137 470.3861 0.9988451
b[18] 1.1830804 0.07722528 1.0549612 1.3046221 397.6879 0.9979990
b[19] 1.9731130 0.07944392 1.8428377 2.0924031 509.8644 1.0024707
b[20] 1.6173833 0.08085352 1.4745344 1.7402543 361.8035 1.0021741
b_mu  1.4483017 0.08322592 1.3281445 1.5826809 541.2091 1.0039075
b_sd  0.3454524 0.06676705 0.2545205 0.4594744 279.2551 0.9981472
sigma 0.2301938 0.04110536 0.1732079 0.3027420 114.1453 0.9980733

a & b vary by clusters, but are correlated

precis(m2, depth=3)
                  mean         sd       5.5%       94.5%     n_eff     Rhat4
b[1]         1.0143803 0.08744104  0.8837770  1.15801625  721.4040 0.9987863
b[2]         1.4569767 0.08179049  1.3351552  1.59187810  536.6146 1.0006891
b[3]         1.2426888 0.08137205  1.1261775  1.37907975  620.6789 0.9994441
b[4]         2.0800591 0.08071669  1.9575744  2.22099455  650.7942 0.9992440
b[5]         1.2136164 0.07634753  1.0985469  1.32535325  629.5236 0.9999293
b[6]         1.2392340 0.07849195  1.1139775  1.35824845  486.3734 0.9988040
b[7]         1.8445091 0.08234411  1.7169503  1.97709040  704.2020 1.0004837
b[8]         1.3839925 0.07384597  1.2634173  1.49657395  577.2652 0.9980004
b[9]         0.8851179 0.07599035  0.7585436  1.00713860  487.2414 0.9980013
b[10]        1.3015516 0.07936911  1.1725432  1.43605375  598.3270 0.9985473
b[11]        1.5967882 0.07801285  1.4753589  1.71759745  527.3257 1.0000229
b[12]        1.3017448 0.07661439  1.1846277  1.42124360  673.8910 0.9984481
b[13]        1.3535050 0.07775585  1.2292371  1.47398155  657.0217 0.9988271
b[14]        2.0786186 0.08421506  1.9410978  2.20668740  472.8595 0.9987464
b[15]        1.3379790 0.07069365  1.2325435  1.45309235  529.5274 0.9980077
b[16]        1.4011446 0.07276654  1.2914682  1.51331075  649.7329 0.9980513
b[17]        1.4734786 0.07084068  1.3650928  1.58876935  750.1340 0.9985369
b[18]        1.1817457 0.07258477  1.0719075  1.29910425  484.5213 0.9980433
b[19]        1.9719898 0.08109897  1.8410768  2.09867885 1118.9537 0.9992206
b[20]        1.6216965 0.07526039  1.5038094  1.73946045  543.1927 0.9990044
a[1]         5.6774484 0.29440941  5.2285598  6.11345725  766.1885 0.9981898
a[2]         4.6985481 0.27857557  4.2427601  5.13432750  601.9208 1.0005713
a[3]         6.2984875 0.26145359  5.8670361  6.69775910  613.5075 0.9989312
a[4]         4.6387087 0.26812261  4.1674365  5.06811355  550.1051 0.9990865
a[5]         7.5493572 0.26083686  7.1196535  7.94289415  692.7991 0.9984522
a[6]         4.6603255 0.26946141  4.2316747  5.09041660  518.9945 0.9996686
a[7]         3.8441048 0.28784602  3.3698946  4.29179275  759.3335 0.9992679
a[8]         2.9212086 0.25510791  2.5156858  3.31248930  577.8133 0.9980125
a[9]         5.6906856 0.26501390  5.2771532  6.10849150  520.9542 0.9981407
a[10]        5.1970856 0.28115695  4.7533679  5.63334530  671.4429 0.9981073
a[11]        5.8557941 0.26654734  5.4458153  6.27611195  551.3730 1.0046470
a[12]        4.1716393 0.26398835  3.7553513  4.57763940  678.5736 0.9980023
a[13]        4.4036094 0.25606265  3.9938437  4.83258585  623.0677 0.9980008
a[14]        1.1732646 0.28256899  0.7603441  1.63746960  459.9166 0.9979982
a[15]        5.2541280 0.24457845  4.8640175  5.66460605  516.2838 1.0005499
a[16]        5.8344220 0.25412395  5.4427746  6.24292575  749.5599 0.9990374
a[17]        4.0016485 0.24891853  3.6205492  4.38283615  802.1613 0.9985174
a[18]        4.2647436 0.26098545  3.8904597  4.68331845  624.8952 0.9985104
a[19]        5.7149885 0.29141917  5.2949689  6.20065335 1091.2689 0.9982974
a[20]        3.5101416 0.26644409  3.0791245  3.95307420  584.6252 0.9983197
a_mu         4.7318570 0.30892371  4.2267149  5.20004140  777.5271 0.9981810
b_mu         1.4486811 0.07625235  1.3275034  1.57009370  640.7886 0.9986301
sigma_ID[1]  1.3803440 0.21540617  1.0699386  1.73998950  389.9947 0.9980070
sigma_ID[2]  0.3493480 0.05968775  0.2621147  0.45054538  483.9477 0.9986399
Rho[1,1]     1.0000000 0.00000000  1.0000000  1.00000000       NaN       NaN
Rho[1,2]    -0.3844072 0.17204899 -0.6324356 -0.07074392  491.9179 0.9988999
Rho[2,1]    -0.3844072 0.17204899 -0.6324356 -0.07074392  491.9179 0.9988999
Rho[2,2]     1.0000000 0.00000000  1.0000000  1.00000000       NaN       NaN
sigma        0.2233401 0.03795678  0.1694872  0.29153988  154.1904 1.0168633

Did we estimate the correlation correctly?

Recovering the parameters

red is estimate without correlation blue is estimate with correlation

A non-centered versions

m3 <- ulam(
  alist(
    V ~ dnorm(mu, sigma), 
    mu <- (a_mu + alpha[ID, 1]) + (b_mu + alpha[ID, 2])*time,
    
    # priors
    # adaptive priors - non-centered
    transpars> matrix[ID, 2]:alpha <-
      compose_noncentered( sigma_ID, L_Rho_ID, z_ID),
    
    matrix[2,ID]:z_ID ~ normal( 0 , 1 ),
    a_mu ~ dnorm(4,1.5),
    b_mu ~ dnorm(1,1),
    vector[2]:sigma_ID ~ dexp(2),
    
    cholesky_factor_corr[2]:L_Rho_ID ~ lkj_corr_cholesky(2),
    sigma ~ dexp(2),
    
    # compute ordinary correlation matrices from Cholesky factors
    gq> matrix[2,2]:Rho_ID <<- Chol_to_Corr(L_Rho_ID)
  ), data = df
)
               mean       sd       5.5%       94.5%    n_eff     Rhat4
Rho[1,1]  1.0000000 0.000000  1.0000000  1.00000000      NaN       NaN
Rho[1,2] -0.3844072 0.172049 -0.6324356 -0.07074392 491.9179 0.9988999
Rho[2,1] -0.3844072 0.172049 -0.6324356 -0.07074392 491.9179 0.9988999
Rho[2,2]  1.0000000 0.000000  1.0000000  1.00000000      NaN       NaN
                  mean        sd       5.5%     94.5%    n_eff    Rhat4
Rho_ID[1,1]  1.0000000 0.0000000  1.0000000 1.0000000      NaN      NaN
Rho_ID[1,2] -0.3795035 0.2037607 -0.6572971 0.0284156 122.0967 1.042638
Rho_ID[2,1] -0.3795035 0.2037607 -0.6572971 0.0284156 122.0967 1.042638
Rho_ID[2,2]  1.0000000 0.0000000  1.0000000 1.0000000      NaN      NaN